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' Abstract. We report a new pulsar population synthesis 
based on Monte Carlo techniques, aiming to estimate the 
contribution of galactic radio pulsars to the continuous 
. gravitational wave emission. Assuming that the rotation 
' periods of pulsars at birth have a Gaussian distribution, 
. we find that the average initial period is 290 ms. The num- 
ber of objects with periods equal to or less than 0.4 s, and 
therefore capable of being detected by an interferometric 
gravitational antenna like VIRGO, is of the order of 5100- 
7800. With integration times lasting between 2 and 3 yr, 
' our simulations suggest that about two detections should 
. be possible, if the mean equatorial ellipticity of the pul- 
sars is £ =10^^. A mean ellipticity an order of magnitude 
. higher increases the expected number of detections to 12- 
' 18, whereas for e < 10~^, no detections are expected 
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1. Introduction 

The first large wide-band gravitational interferometric 
antennas, the French-Italian VIRGO and the American 
LIGO, should be fully operational at the beginning of the 
next century. The expected sensitivity of these detectors in 
terms of the gravitational strain is of the order of hw 10"^^ 
for short-lived, impulsive bursts, and h« 10^^^ for peri- 
odic "long-lived" sources for which integration times of 
about 10^ s are possible. Clearly, a major concern of these 
detectors is to understand and minimize the sources of 
noise which define the sensitivity of the antenna, but the 
identification of a possible "standard" source (or sources), 
having a well defined waveform, is also a relevant problem 
in the context of the signal detection. 

Among the possible candidates, it is expected that neu- 
tron stars, which are known to be quite abundant in the 
Galaxy, could rank aamong conspicuous emitters of grav- 
itational waves (GW). Rotating neutron stars may have a 
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time-varying quadrupole moment and hence radiate GW, 
by either having a triaxial shape or a misalignment be- 
tween the symmetry and the spin axes, which produces a 
wobble in the stellar motion (Ferrarife Ruffini 1969; Zim- 
merman & Szedenits 1979; Shapiro & Teukolsky 1983; 
de Araujo et al. 1994). Moreover, fast rotating proto- 
neutron stars may develop different instabilities such as 
the so-called Chandrasekhar-Friedman - Schutz (CFS) in- 
stability (Chandrasekhar 1970; Friedman & Schutz 1978), 
responsible for the excitation of density waves traveling 
around the star in the sense opposite to its rotation, or 
to undergo a transition from axi-symmetric to triaxial 
shapes through the dynamical "bar" -mode instability (Lai 
& Shapiro 1995). All these mechanisms are potentially 
able to emit large amounts of energy in the form of GW. 

Continuous GW emission by the entire population 
of rotating neutron stars. Assumed to have small devi- 
ations from axisymmetry, raises the possibility of detect- 
ing their combined signals. The contribution of individual 
observed sources, like galactic pulsars, has been consid- 
ered by Barone et al. (1988), Velloso et al. (1996), among 
others. However, detected radio pulsars constitute only a 
small fraction of the actual galactic population. Attempts 
to estimate the gravitational strain of the total pulsar 
population have been made by Schutz (1991), Giazotto, 
Bonazzola & Gourgoulhon (1997, hereafter GBG97), de 
Freitas Pacheco & Horvath (1997), Giamperi (1998). Ex- 
cepting Schutz (1991), all these authors proposed to de- 
tect the square of the gravitational strain, and the side- 
real modulation of the integrated signal was examined by 
GBG97 as well by Giamperi (1998). In order to perform 
such an estimation, the actual distribution of the rotation 
periods must be known. The rotation period is a critical 
parameter, because the gravitational strain depends on 
the square of the angular velocity. Unfortunately, the ob- 
served distribution does not necessarily reflect the real pe- 
riod distribution due to different selection effects present 
in all pulsar searches. GBG97 assumed for their estimate 
an (optimistic) average rotation period of 5 ms, which 
now seems a rather short value. If we exclude binary mil- 
lisecond pulsars, which have probably been spun-up by 
accretion mechanisms, several observational facts indicate 
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that population I pulsars are born with smaller rotation 
velocities (Bhattacharya 1990). Moreover, several popula- 
tion synthesis calculations also suggest higher initial peri- 
ods (Narayan 1987; Bhattacharya ct al. 1992). Magnetic 
coupling between the proto-neutron star and the outer en- 
velope could be an efficient mechanism for transferring an- 
gular momentum, so reducing the initial angular velocity. 
Proto-neutron stars are expected to be hot, with temper- 
atures around 10® K. In this situation, depending on the 
viscosity, a rapidly rotating proto-neutron star may excite 
r-modes, emitting GW which decelerates considerably the 
star. According to computations by Andersson, Kokkotas 
& Schutz (1999), this mechanism sets a limit around P « 
20 ms for the fastest newly born neutron stars. However, 
it should be emphasized that this value is a function of 
the adopted (uncertain) viscosity law. 

In the present work we revisit the problem of the con- 
tribution to the continuous GW emission from the pulsar 
population inside the galactic disk. We estimate the actual 
period distribution by using population synthesis based 
on Monte Carlo methods, looking for the best pulsar pop- 
ulation parameters by comparing our simulation results 
with data. Using the model population, we compute the 
gravitational strain for each pulsar emitting GW in the 
frequency band of VIRGO, the total square of the signal 
as well the amplitude modulation due to the sidereal mo- 
tion. We show, in agreement with the expectations made 
by dc Freitas Pachcco & Horvath (1997), that the signal 
is strongly dominated by a few fast pulsars and/or by the 
nearest ones. The plan of this paper is the following: in 
section 2 we describe our population synthesis method; in 
section 3 we present the properties of the simulated pop- 
ulation; in section 4 we estimate the contribution to the 
continuous gravitational strain of such a population and 
finally, in section 5 we present our main conclusions. 

2. The Population Synthesis 

2.1. The Rotation Period Evolution 

The evolution of the rotation period of a neutron star, in 
the framework of the magnetic dipole model, is given by 
(Pacini 1968) 

where /i = ^BR^ is the surface magnetic moment, B is 
the magnetic field, Rp and I are respectively the radius 
and the moment of inertia of the neutron star, c is the 
velocity of light and a is the angle between the rotation 
axis and the magnetic dipole. Integration of this equation 
with constant a gives the rotation period evolution of the 
so-called standard model. 

In the absence of abrupt changes of the magnetic 
torque or of the internal structure of the star, solutions 
of eq. (1) indicates a smooth increase of the rotation pe- 
riod with time. However, most pulsars do not slow down 



regiilarly, but display variations in their spin rates in the 
form of glitches (Shemar & Lyne 1996) and timing noise. 
An important aspect of glitches is the persistent increase 
of the spin-down rate following these events, which could 
be due to a sudden and permanent increase of the exter- 
nal torque (Link & Epstein 1997) .In order to explain the 
glitches observed in the spin evolution of the Crab pulsar, 
Alpar & Pines (1993) suggested a reduction of the moment 
of inertia, induced by changes in the internal structure of 
the star. This would produce a spin-up of the star due 
to angular momentum conservation, but a decrease in the 
spin-down rate, contrary to observations. As the star spins 
down, it becomes less oblate, inducing stresses that may 
lead to starquakes if the yield strength of the solid crust is 
exceeded. Starquakes may affect the braking mechanism if 
they are able to change the position of the magnetic mo- 
ment with respect to the rotation axis. If the structural 
relaxation after a starquake occurs asymmetrically about 
the rotation axis due, for instance, to magnetic stresses, 
the figure and spin axes may become misaligned. Under 
this condition, the star precesses and relaxes to a new 
equilibrium state, corresponding to a new orientation of 
the magnetic dipole with respect to the rotation axis. This 
possibility was examined in a recent work by Link, Franco 
& Epstein (1998). In their picture, a spin-down pulsar 
reduces its equatorial radius by shearing material across 
the equator and moving material along faults to higher 
latitudes. Such crustal motions may produce an increas- 
ing misalignment between the rotation and magnetic axes, 
providing a natural explanation for the observed increases 
in the spin-down rate following glitches in the Crab, PSR 
1830-08 and PSR 0355-1-54. 

The migration of the magnetic dipole fromm an initial 
angle a to an orthogonal position with respect to the spin 
axis was also considered by Allen & Horvath (1997). They 
showed that with such a variable torque, non-canonical 
braking indices may be explained. Here we adopt the sce- 
nario developed by Link, Franco & Epstein (1998), assum- 
ing that the magnetic dipole migrates from an arbitrary 
initial angle to a final orthogonal position in a timescale 
ta- For our present purposes, this migration is modeled by 
the equation 

sm2(a(t)) = 1 - noe"*/*" (2) 

We emphasize that this equation is not the result of any 
physical model, but only a phenomenological representa- 
tion of the above picture, in the sense that a given pulsar 
born with an initial angle ao = acos{^Jno), will reach an 
orthogonal configuration in a timescale ta ■ 

Combining eq. (1) and (2) and integrating, we obtain 
for the evolution of the rotation period 

P = Po[l + l-no^(l-e-*/*")]V2 (3) 

To To 

where Pq is the initial rotation period and Tq = ^2 j^e 
is the magnetic braking timescale. For t >> ^q, we re- 
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cover essentially the standard model evolution. However, 
eq. (3) allows one to explain the existence of objects with 
short periods and with a relatively small deceleration rate, 
which would be the consequence of a small initial misalign- 
ment between the magnetic and the spin axes, but not due 
to a small magnetic field strength. Note that in our pic- 
ture these pulsars are young, in spite of having a large 
indicative age = 

Another important parameter characterizing the evo- 
lution of the rotation period is the so-called braking index 
N defined as 



For the standard magnetic dipole model, the expected 
braking index is TV = 3. However, the few pulsars with N 
measured accurately show deviations from the Standard 
model; 7V=2.518 for the Crab, 7V=2.837 for PSR1509-58, 
N^2M for PSR0540-69. The exceptionally low value for 
the Vela pulsar N = 1.40 was recently reported, suggesting 
that the braking torque is not due to a "pure" stationary 
dipole field (dc Frcitas Pachcco & Horvath 1997). From 
eqs. (3) and (4), the predicted braking index for our model 
is 

N=[3-2{*f)ctg\am (5) 

This equation indicates that for young pulsars with ages 
less than ta the braking index N may be less than three, 
tending to the canonical value when t » ta- 

The angle a may vary " discontinuously" when the 
crust suddenly cracks or almost "continuously" as eq. (2) 
suggest. It is worth mentioning that pulsars like PSR1509- 
58 and PSR0540-69 have not yet displayed any glitch ac- 
tivity even though their braking indices are less than the 
canonical value N = 3. This could be an indication that 
other braking mechanisms are operative such as, for in- 
stance, anomalous electro-magnetic torques recently dis- 
cussed by Casini & Montemayor (1998). 

For the Crab pulsar, since the true age, the braking in- 
dex N and are known, we can compute using the equa- 
tions above the dipole migration rate < ctga-^ >, which 
is about 9.2x 10~^ rad/yr for a solution where the present 
angle between the spin and the dipole axes is about 64°. 
Prom glitch data spanning a time interval of about 20 yr, 
Allen & Horvath (1997) find a migration rate 1.8x10""'^ 
rad/yr, which is about five times smaller than the pre- 
dicted value. However, as those authors remarked, such a 
time interval is still very short and discrepancies should 
be expected when a comparison is performed between the 
migration rate derived from the average of few events and 
that obtained from a continuous variation model. For the 
other remaining three pulsars, the average "continuous" 
migration rate is < ctga ■ ^ > = 6.5x 10~^ rad/yr (Allen 
& Horvath 1997). These rates are comparable to within 
an order of magnitude, reinforcing the idea of a common 



origin and giving some support to our adopted spin-down 
model. 

2.2. The Method 

Initial spin periods and magnetic fields are important pa- 
rameters tracing the formation history of neutron stars 
and their interaction with the surrounding medium. The 
pulsar population study by Narayan (1987) raised the so- 
called "injection" problem, which would connect the initial 
period and magnetic field in the sense that slow born neu- 
tron stars would be associated with high field strengths. 
Unfortunately, the ages of most pulsars are unknown and 
present measurements for estimating the period and its 
time derivative are insufficient in order to estimate the 
initial field and spin period. In this case, methods of pop- 
ulation synthesis can be useful tools to estimate statis- 
tical values of these parameters, assuming ad-hoc distri- 
bution laws. Moreover, the observed pulsar population is 
strongly dominated by luminous objects, and since the 
radio emission depends on the spin period as well as on 
the spin down rate, present samples are biased, exclud- 
ing long-period and high-field objects as well as old and 
faint pulsars. Comparison with actual data requires the 
inclusion of these selection and observational bias in the 
generation of synthetic catalogs, which are more treatable 
when Monte Carlo techniques are employed. 

In order to generate a synthetic sample of pulsars, we 
need to model: (a) their Distribution in the Galaxy, (b) the 
initial parameters aflfecting their evolution, (c) the effects 
of the interstellar medium in their detection and, (d) the 
sensitivity of a given survey. 

The distribution in the Galaxy was modeled by as- 
suming a constant birthrate until a maximum age tmax, 
considered as a free parameter in our numerical experi- 
ments. The ratio between the total number Np of ggener- 
atc pulsars and tmax gives the mean birthrate. For a given 
experiment, Np is the number of input pulsars required to 
reproduce the number of objects in the data sample, when 
observational and selection constraints arc imposed. 

Pulsars are strongly concentrated in the galactic plane 
and their distribution along the z-axis can be represented 
by an exponential law with a scale of height of about 100 
pc. Some previous studies suggest a possible correlation 
between the height above the galactic plane and age or 
magnetic field (Narayan & Ostriker 1990), but here we 
assume that the initial z-coordinate is a statistically in- 
dependent variable. In the galactic plane, the density of 
pulsars probably decays exponentially, in the same way 
as population I objects. Recent galactic mass profile mod- 
els, based on infrared data, suggest a short distance scale 
of about 2.3 kpc (Ruphy et al. 1996); this was adopted 
in our numerical experiments. In order to compute helio- 
centric distances, the galactocentric distance of the sun 
was assumed to be 8.5 kpc. Some pulsars have significant 
proper motions, implying large initial velocities that must 
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be taken into account when computing distances, since 
they affect tfie orbital motion in the Galaxy. An isotropic 
Gaussian velocity distribution with 1-D velocity disper- 
sion of 100 km/s (Lorimcr ct al. 1993) was adopted to 
simulate the pulsar kick velocity, but our results are not 
significantly changed if the velocity dispersion is increased 
up to 200 km/s. Once the initial R and z coordinates for a 
given pulsar were established, if its kick velocity is zero, we 
assumed that the pulsar follows a circular orbit using the 
galactic rotation curve by Schmidt (1985). For non-zero 
kick velocity, the new energy and angular momentum of 
the orbit were computed and the object was settled in a 
new trajectory. The heliocentric distance is computed af- 
ter a time t (the age of the pulsar), taking into account 
the differential rotation between the sun and the pulsar. 
The present galactic coordinates ( 1, b) of the object arc 
also computed, since these are required to establish if the 
pulsar should be included or not in the area of the sky 
covered by a given survey, as we shall see later. 

In the present simulations we have assumed that the 
initial rotation period of pulsars can be represented by a 
Gaussian distribution, characterized by a mean value Po 
and a dispersion op^. We have assumed that the break-up 
velocity is about 12000 rad/s, thus pulsars generated with 
rotation periods less than 0.5 ms were not included in the 
simulated sample. These parameters were allowed to vary 
in order to verify their effects on the results of different 
numerical experiments. A similar procedure was adopted 
for the parameter tq (see, eq. 3), which is related to the 
magnetic field of the pulsar, for which a log-normal dis- 
tribution was assumed. Simulations by Bhattacharya et 
al. (1992) and by Mukherjee & Kembhavi (1997) suggest 
that the timescale for magnetic field decay is longer then 
the pulsar lifetime {td > 100-160 Myr), thus we have as- 
sumed a constant field in our numerical experiments. The 
initial angle ao (or, equivalently, the value of no) between 
the magnetic dipole and the spin axes is calculated sup- 
posing a random distribution. Once these parameters are 
assigned for a given pulsar {ta is kept fixed during a giv- 
ing experiment), the period and the time derivative can 
be computed from eq.(3). 

The comparison between a simulated sample with ac- 
tual data requires the evaluation of the expected radio flux 
density S^ at a given frequency (400 Mhz in general). To 
compute Su, one needs to know the distance and the lu- 
minosity of the pulsar. A theory able to predict the pulsar 
radio luminosity is still missing, but the spin period and 
its time derivative are expected to be the main variables 
on which the radio power depends, since the magnetic 
field can be expressed in terms of these variables in the 
standard model. Different power laws have been suggested 
in the literature connecting P and P (Lyne, Ritchings & 
Smith 1975; Stollman 1987; Emmering & Chevalier 1989). 
When these laws are compared with observational data a 
large dispersion is foimd, which can be explained by a bad 
understanding of the radio emission mechanism and/or by 



errors in the distance, estimates of which are derived by 
using the dispersion of radio waves through the interstel- 
lar plasma and by adopting a model for the electron den- 
sity distribution inside the galactic disk. The effects of the 
adopted luminosity law in numerical simulations were dis- 
cussed by Lorimer et al. (1993), who concluded that, at 
the present time there appears to be no satisfactory model 
to account for observed pulsar luminosities. 

The radio luminosity ought to be a decreasing function 
of the spin period, since slow pulsars missing in different 
surveys are probably too weak to be detected. Models as- 
suming a radio luminosity proportional to the cube root of 
the spin-down rate lead to an inverse dependence on the 
period (Narayan 1987), while other proposed power laws 
indicate even higher exponents. In the present work, we 
have adopted a different empirical luminosity law, impos- 
ing the following exponential fate on the radio emission 

Lr = AP°'e-^^ (6) 

The parameters in this relation were determined by an 
optimised fitting procedure, using the upgraded pulsar 
catalog of Taylor, Manchester & Lyne (1993). If the lu- 
minosity is given in mJy.kpc^, and the period is in sec- 
onds, the resulting numerical values of the parameters 
are: A=3.3xl0^, a — 0.34 and 7 = 1.02. These figures are 
slightly different from those given by de Freitas Pacheco & 
Horvath (1997) because here we have used a larger sample 
in the fitting procedure. 

The radio emission is not isotropic, implying that pul- 
sars are not visible from all directions, and this imposes 
another constraint on their detectability. Following the ap- 
proach used by Emmering & Chevalier (1989), which will 
be adopted here, the beaming fraction of the sky covered 
by the pulsar emission cone is 

/ = (l-cos^) + (|-^)sin^ (7) 

where 20 is the aperture of the emission cone. Since short 
period pulsars seem to have wider emission cones, we have 
adopted the beaming model proposed by Biggs (1990), 
which gives for the cone aperture 



where is in degrees and P is in seconds. It is worth 

mentioning that the study by Rankin (1990) of a large 
sample of pulsars also indicates a variation of 9 inversely 
proportional to the square root of the period. 

Our data are a culled sample of 389 pulsars extracted 
from the catalog of Taylor, Manchester & Lyne (1993), 
where only objects that have originated from population I 
stars and with all required parameters measured were con- 
sidered. Pulsars not included in any of the original surveys 
were also discarded. Since this catalog is a compilation of 
different catalogs covering specific sky areas and having 
a well defined sensitivity, the objects were distributed in 



5 



sub-classes, according to the original surveys in which they 
were first detected. 

Each simulated catalog to be compared with these 
data is prepared with the following prescriptions. For each 
pulsar generated with a given set of initial parameters, 
we have followed its evolution in order to compute the 
present P, P values and galactic coordinates, according 
to the scheme described above. Then the radio luminos- 
ity and the flux density S^, = ^ are calculated, as well 
as the probability for the pulsar beam to cross the earth. 
The galactic coordinates define the survey (or surveys) 
included in the general catalog that covered this region. 
Then, we compare the expected flux density with the min- 
imum detectable flux density of that survey through the 
relation (Dewey et al. 1984; Manchester et al. 1996) 



Smin — ^0 {Tj. + Tsfcj/ ) [ 



W 



{P-W)' 



11/2 



(9) 



where W is the observed (broadened) pulse width, T^. and 
Tsky are respectively the system and the sky background 
temperature in the considered line of sight. Aq is a nor- 
malization constant, defined as 



^0 



{S/N)f3 



Gas/ijlpAutint) 



(10) 



In this equation, Ga is the antenna gain, np is the num- 
ber of polarizations used, Ai^ is the receiver bandwidth, 
tint is the integration time, S/N is the signal-to-noise ratio 
and /3 is a numerical factor, including other loss processes. 
The variation of the sky temperature with the galactic 
coordinates was estimated using the relation by Narayan 
(1987). The intrinsic pulse width W,. was assumed to be 
equal to 4% of the period ( We=0.04P) and to model 
the pulse broadening due to instrumental effects, disper- 
sion and scattering, we follow the approach by Dewey et 
al.(1984) (see also Lorimer et al. 1993), writing for the re- 
lation between the observed and the intrinsic pulse width 



(11) 



where Tsam is the effective sampling time, T^ca is the 
broadening due to diffusion through the interstellar 
medium and tum is the channel broadening. We have used 
the relation given by Bhattacharya et al. (1992) for Tsca-, 
namely, 

Tscaims) = 10-4-62+1.14;os(DM) ^ ^^-•A.11^i.mog{DM) (jj) 

where DM is the electron dispersion measure in the line of 
sight of the pulsar in pc.cm"^. For the channel broadening, 

we have used the relation 



tdm{iJ'S) = Cdm-DM 

where, following Stokes et al. (1986) 

BchiMHz) 



C. 



DM 



8.3- 



v^{GHz) 



(13) 
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with Bc/i being the channel bandwidth and v the observing 
frequency. Table 1 gives the main parameters characteriz- 
ing the different catalogs adopted in our simulations. Note 
that values inside braces correspond to high latitude data. 

Our procedure insures that the simulated catalog will 
have the objects distributed in the same proportion as in 
the global catalog, taking into account the specific sensibil- 
ity and the sky area covered by each survey. The numerical 
experiments generate objects according to these prescrip- 
tions until the simulated number of detectable pulsars is 
equal to the number in our data sample. For a given exper- 
iment, defined by a set of initial parameters, the number 
of runs is comparable to the number of objects in the sam- 
ple and the final result is a suitable average, in order to 
avoid statistical fluctuations. 

For each experiment so defined, we have compared 
the resulting distributions with those derived from actual 
data. We searched for optimize the various input param- 
eters that describe the pulsar population, controlling the 
fit quality through and Kolmogorov-Smirnov tests. 

3. Properties of the Simulated Pulsar Population 

The possibility that the magnetic axis could migrate gives 
not only a possible explanation for the non-canonical brak- 
ing index observed in some pulsars, but also a mechanism 
affecting the P and age distributions. In general, objects 
with small deceleration rates (P w lO"'^^) are old and 
have low magnetic fields in the standard picture. In our 
scenario the magnetic torque controlling the deceleration 
rate, depends not only on the field strength but also on the 
angle between the dipole and the spin axes. Therefore, our 
simulated catalog contains some young objects born with 
a small misalignment between those axes, having a high 
field strength but a small deceleration rate. Unfortunately, 
it is not possible to disentangle easily both situations and, 
as a consequence, the derived value of the dipole migration 
timescale ta has a large uncertainty: ta = 10000±4000 yr. 
These values imply migration rates spanning the interval 
(2 - 18)xl0~^ rad/yr, which compare with rates derived 
by other authors (see section 2.1). 

The optimized parameters characterizing the distribu- 
tion of the initial period and that of the magnetic brak- 
ing timescale tq are given in table 2 for tmax equal to 24 
Myr and 100 Myr. These parameters represent the best 
compromise able to fit adequately to data, the simulated 
distributions of period and its derivative, as well as the 
distribution of distances to the sun. Note that the param- 
eters characterizing the Gaussian representing the initial 
distribution of rotation periods are not sensitive to the 
adopted value for tmax, but this is not the case for the 
parameter tq. 

In figure 1 we have plotted in the plane XY, which 
coincides with the galactic plane (the origin of the coordi- 
nates is placed at the galactic center), the observed pulsar 
population and our simulated population. The superpo- 
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Table 1. Survey Parameters 



Survey 


(K) 


Tsam (ms) 


Ci3M(ms.pc ^.cm^) 


Ao(mJy/K) 


Arecibo 1 


110 


33 


0.0094 


0.038 


Arccibo 2 


90 


0.6 


0.0063 


0.085 


Arecibo 3-4 


62 


1.0(0.5) 


0.0082(0.026) 


0.02(0.032) 


Green Bank 1 


170 


33 


0.28 


0.30 


Green Bank 2 


30 


33 


0.28 


0.13 


Green Bank 3 


30 


4.0 


0.035 


0.18 


JodrcU Bank 


110 


80 


0.30 


0.36 


Molonglo 2 


210 


40 


0.06(0.3) 


0.17 


Parkes 


30 


0.60 


0.0125 


0.19 



sition of both populations is quite good, suggesting that 
selection effects were adequately taken into account by the 
model. In figure 2 wc compare the observed period distri- 
bution with simulated data, where error bars correspond 
to the dispersion derived from five hundred experiments. 
Figure 3 shows the comparison between observed and sim- 
ulated data for the period derivative, whereas figure 4 dis- 
play both heliocentric distance distributions. These figures 
correspond to simulations with tmax = 24 Myr and they 
illustrate the fit quality. 

The last line of table 2 gives the average magnetic 
field of the simulated observed population, which is about 
2.5x10^^ G, in agreement with values derived from a di- 
rect application of the "standard" pulsar model. However, 
the average field of the true or " unseen" population is one 
order of magnitude higher, namely, 2.5x10^^ G. Most of 
these high-field pulsar have rather long periods and thus 
are radio-quiet. The relevance of this high- field pulsar pop- 
ulation in the context of magnetars will be discussed in a 
forthcoming paper (Regimbau & de Freitas Pacheco 2000). 

Our simulations indicate that the data are better ex- 
plained if an initial distribution of periods is assumed in- 
stead of an unique (mean) initial value. Bhattacharya et 
al. (1992) assumed an initial period equal to 100 ms, com- 
parable with the average value derived from our simula- 
tions (290 ms). Models by Lorimer et al. (1993) are not 
directly comparable, since they have assumed a magnetic 
field decay with a timcscalc of 10 Myr. For the present 
purposes, we will focus our attention on the pulsar pop- 
ulation with periods less than 0.4 s, which contributes to 
the continuum GW emission from the galactic disc. This 
upper bound is expected to be attained in a second phase 
of the VIRGO experiment. In spite of the bulk of the pop- 
ulation having higher periods, the dispersion in The ini- 
tial periods guarantees the existence of objects in that 
frequency band. Prom our simulations, we predict about 
60-90 (single) pulsars with P < 80 ms in the Galaxy, a 
number not in contradiction with present data. The pe- 
riod distribution for pulsars satisfying the condition P < 
0.4 s is shown in figure 5, and their heliocentric distance 
distribution is given in figure 6. The number of this sub- 



population is in the range 5100 - 7800, according to our 
simulations performed with imax equal to 24 Myr and 100 

Myr respectively. The contribution of this population to 
the gravitational strain is estimated in the next section. 

4. The Gravitational Strain 

4.-1- The Equations 

As we have mentioned above, pulsars could emit GW by 

having a time-varying quadrupole moment produced ci- 
ther by a slight asymmetry in the equatorial plane ( as- 
sumed to be orthogonal to the spin axis) or by a mis- 
alignment between the symmetry and angular momentum 
axes, case in which a wobble is induced in the star mo- 
tion. In the former situation the GW frequency is equal 
to twice the rotation frequency, whereas in the latter two 
modes are possible: one in which the GW have the same 
frequency as the rotation, and another in which the GW 
have twice the rotation frequency. The first mode domi- 
nates by far at small wobble angles while the importance 
of the second increases for large values. 

Here wc neglect the possible preccssional motion and, 
in this case, the two polarization components of GW emit- 
ted by a rotating neutron star are (Zimmerman & Szeden- 
its 1979; Bonazzola & Gourgoulhon 1996 ) 

h+{t) = 2A{1 + cos^ i) cos(20t) (15) 

and 

h^{t) = AAcosi.sin{2Q.t) (16) 

where i is the angle between the spin axis and the wave 
propagation vector, assumed to coincide with the line of 
sight, 

A = -^eh^n^ (17) 

G is the gravitation constant, c is the velocity of light, r 
is the distance to The source, O is the angular rotation 
velocity of the pulsar and the ellipticity e is defined as 
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Table 2. Pulsar Population Properties 



Parameters 


tmax = 24 Myr 


tmax = 100 Myr 


Po (ms) 


290±20 


290±20 


CTpj, (ms) 


90±20 


90±20 


In(ro) (yr) 


9.0±0.5 


11.0±0.5 


ffin(To) (yr) 


3.6±0.2 


3.4±0.2 


< logB > (G) 


12.4 


12.1 



with the lij being the principal moments of inertia of the 
star. 

The detected signal by an interferometric antenna is 

h{t) = h+{t)F+ {e, cp,^)+hy^ {t)Fy (6, 4>, V) (19) 

where and Fx are the beam factors of the interfer- 
ometer, which are functions of the zenith distance 9, the 
azimiith (f> as well as of the wave polarization plane orien- 
tation ip. Notice that the angles 6 and are functions of 
time due to the Earth's rotation, introducing a modula- 
tion of the signal. The explicit functions for the beam fac- 
tors, taking into account the geographic localization and 
the orientation of the VIRGO antenna were taken from 
Jaranowski, Krolak & Schutz (1998). 

Concerning the detection strategy, recall that the pop- 
ulation derived from our simulations of potential GW 
emitters is in the range 5100 - 7800. In this case, according 
to the conclusions by GBG97, it becomes advantageous to 
search for individual detections instead of the total square 
amplitude. Here we have simulated both strategies. In the 
former case, the strain amplitude was calculated for each 
pulsar satisfying the condition P < 0.4s, using eq.(19) and 
assuming a random orientation for the inclination i of the 
spin axis as well for the orientation of the polarization 
angle ip. Since the signal is modulated at twice the rota- 
tion frequency, in general much shorter than the sidereal 
period, we have assumed also a random phase for the rel- 
ative orientation of the detector with respect to the equa- 
torial coordinate system. In the other case, the procedure 
adopted to calculate the total square of the strain was the 
following: first, we have squared equ.(19) and then per- 
formed an average in the time interval r, satisfying the 
condition, 0.4 s < r < one day. In this case, the cross 
terms give a null contribution and we obtain 

< h'^it) >= A^[2{1 + cos^ if Flit) + 8 cos^ iF^ (t)] (20) 

This equation was used to compute the contribution 

from all simulated objects satisfying P < 0.4 s, assuming 
again a random orientation for i and ip. 

4.2. The Results 

The main differences between the present approach and 

previous calculations should be emphasized. Our proce- 
dure allows a more realistic estimate of the rotation period 



distribution, as well the number of pulsars able to con- 
tribute to the gravitational strain. Additionally the spatial 
distribution of those pulsars throughout the galactic disc 
changes for each simulation, although their average prop- 
erties remaining constant. It is thus preferable to present 
the statistics of our numerical experiments, from which 
is possible to estimate the probability of having a signal 
above a given threshold. 

In figure 7 we give the statistics for the gravitational 
strain h. For each experiment we have computed the dis- 
tribution of values of h and then averaged the results of 
500 experiments. Error bars indicate the rmsd for each bin 
of thickness Alog{\ h |) = 0.20. The resulting distribution 
can be represented quite well by a Laplace law, namely 

L{x) = e(-2|-+29-7|) (21) 

where x = log(| h \). This function corresponds to calcu- 
lations performed with an ellipticity e = 10^^. For other 
values of e, it is enough to replace the constant in eq.(21) 
by 23.7 log(£). The probability per pulsar tohave a signal 
above xq is 

poo 

p{xo) = / L{x)dx (22) 

For a continuous source, VIRGO (or LIGO) will be able to 
detect amplitudes of the order of 10~^^, with integra- 
tion times of about 2-3 yr. From the equation above and 
the population number estimated previously, one should 
expect to detect about 2 objects if e = 10~^ or about 12- 
18 pulsars if e = 10~^. If the average ellipticity is smaller 
than 10~^, most of the objects will be below the detection 
threshold, at least for the present antenna sensibility. 

In figure 8 we present the statistics for the square of 
the amplitude, resulting from the different spatial distri- 
bution of the pulsars in each numerical experiment. We 
emphasize that the amplitude of the signal is due essen- 
tially to a few pulsars, in agreement with the conclusions 
derived from the statistics of single objects and with the 
analytical study by de Freitas Pacheco & Horvath (1997). 
As a consequence, the sidereal modulation is not fixed 
by the galactic center-anticenter asymmetry, but by those 
few dominant objects. No typical modulation curve was 
obtained, since the relative positions of these pulsars vary 
from experiment to experiment. Using the equations by 
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Regimbau.: GW from Pulsars 



GBG97 for the signal-to-noise ratio, one should expect 
to detect a signal of < >« 10"'*^ with the presently 
planned VIRGO sensibility. Our simulations indicate that, 
a e = 10^^ a signal of such an amplitude has weak prob- 
ability to be detected (about 1/100) and signals of the 
required amplitude can only be obtained if the average 
pulsar ellipticity is of the order of 10"'*. 

5. Conclusions 

Estimates of the GW emission from the total pulsar pop- 
ulation in the Galaxy require knowledge of the period 
and distance distributions. Selection and detection bias 
are strong and population synthesis methods are neces- 
sary in order to model these effects, to recover the "true" 
population properties from the observed distributions. 

Pulsars are probably born with a range of periods and 
magnetic fields. Here we have assumed ad-hoc initial distri- 
butions for these parameters, which were optimized in or- 
der to reproduce actual data. Our numerical experiments 
suggest that the mean initial period at the pulsar birth 
is equal to 290 ms, confirming the conclusions of other 
simulations. Nevertheless, we predict about 60-90 single 
pulsars in the galaxy with periods less than 80 ms, as a 
consequence of the dispersion in the initial periods. 

The pulsar population satisfying the condition P<0.4 
s, able to contribute to the GW emission within the pass- 
band of VIRGO, amounts to about 5100-7800 objects. 
These numbers are considerably smaller than the estimate 
adopted by GBG97 (1.4x10^). Those authors assumed 
that the population with P< 0.4 corresponds to about 
28% of the total population. The adopted fraction is that 
derived from cataloged pulsars. However, our simulations 
indicate that this fraction is only about 3.5%, exemplify- 
ing the consequences of the detection bias present in all 
surveys. 

VIRGO will be able to detect a gravitational strain 
amplitude of h« 10~^^ with three years of integration. 
Under these conditions, our simulations indicate 2 detec- 
tions if e = 10~^ and up to 12-18 detections if e = 10~^. 
If the average ellipticity is smaller than 10~^, no detec- 
tions are expected, at least with the presently planned 
antenna sensitivity. The total square amplitude resulting 
from this population is below the VIRGO threshold and 
a detectable signal would be produced only if the aver- 
age ellipticity is at least of order of 10"^, a value which 
seems to be excluded by a recent re-analysis of the mag- 
netic and gravitational torques in some (observed) young 
pulsars (Palomba 1999). 

We emphasize that the present results concern 
uniquely the radio pulsar population and the contribu- 
tion of a binary millisecond population to the continuous 
GW emission remains to be estimated. These simulations 
are more difficult, since the rotation period evolution has 
a complex history, including decelerating and accelerating 
phases, which will be discussed in a future paper. 
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Figure Captions 
Figure 1 

Spatial distribution in the galactic plane (origin at the 
galactic center) of both observed (stars) and simulated 
(circles) populations of radio pulsars. The position of the 
sun is also indicated in the diagram. 
Figure 2 

Distribution of rotation periods. Bars represent observed 

binned data and filled circles indicate the result of our 
simulations. Error bars give the rmsd after averaging 500 
numerical experiments. 
Figure 3 

Distribution of period derivatives. Symbols have the same 
meaning as in figure 1. 
Figure 4 

Distribution of hehocentric distances. Symbols have the 
same meaning as in figure 1. 
Figure 5 

Distribution of periods for all pulsars satisfying P<0.4 s. 
Figure 6 

Distribution of heliocentric distances for the pulsar popu- 
lation satisfying P<0.4 s. 
Figure 7 

Distribution of the gravitational strain h for e = 10~^. 
Error bars indicate rmsd after averaging 500 numerical 
experiments. 
Figure 8 

Distribution of the total h^ for e = 10~^, corresponding 
to 7000 numerical experiments. 
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